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ABSTRACT 

Aims. The early acceleration of stellar winds in massive stars is poorly constrained. The scattering of hard X-ray photons emitted 
by the pulsar in the high-mass X-ray binary Vela X-l can be used to probe the stellar wind velocity and density profile close to the 
surface of its supergiant companion HD 77581. 

Methods. We built a high signal-to-noise and high resolution hard X-ray lightcurve of Vela X-l measured by Swift/BAT over 300 
orbital periods of the system and compared it with the predictions of a grid of hydrodynamic simulations. 

Results. We obtain a very good agreement between observations and simulations for a narrow set of parameters, implying that the 
wind velocity close to the stellar surface is twice larger than usually assumed with the standard beta law. Locally a velocity gradient of 
P ~ 0.5 is favoured. Even if still incomplete, hydrodynamic simulations are successfully reproducing several observational properties 
of Vela X-1. 

Key words. X-rays: Binaries, Hydrodynamics, Stars: winds, Accretion, Stars: individual: Vela X-l 


1. Introduction 


The density and velocity of stellar winds in massive stars are 
poorly constrained close to the stellar surface. Eclipsing high- 
mass X-ray binaries with short orbital orbits allow to probe the 
conditions of the inner stellar wind in situ. Thanks to its high X- 
ray flux. Vela X-l is a perfect laboratory, with a pulsar accreting 
and photo-ionizing the wind of its massive (B 0.5 lb) companion. 
The X-ray emission produced close to the neutron star traces 
instabilities of the accretion flow on large scales, which can be 
studied thanks to the variability of the X-ray flux and absorbing 
column density (Walter & Zurita Heras [2007| >. 

The neutron star of Vela X-l weights M, ViV = 1.86 Mg 
( |Quaintrell et al.|[2003[ ) and orbits its massive companion with 

a period of 8.96 days in an almost circular orbit (a =1.76 R t , 
e ~ 0.09; |Bildsten et al.|1997[ >. The wind from its massive stellar 
companion is characterized by a mass-loss rate of ~ 4 x 1CL 6 
M g yr~' (Nagase et ah 1986) and a wind terminal velocity of 
v.„ ~ 1700 km s~* ( Dupree et al.||l980 i, translating to a local 
wind velocity of 400 km/s at 1.2 stellar radius, assuming the 
commonly used [j = 0.5 velocity gradient. The typical X-ray 
luminosity is of the order of ~ 4 x 10 36 erg s _1 , although high 
variability can be observed (Kreykenbohm et al.|1999). 


In this paper we present the hard X-ray lightcurve of Vela X- 

1, folded along the orbit, measured with an unprecedented signal 

to noise thanks to almost eight years of Swift/BAT observations, 

and compare it with the predictions of 2-D hydrodynamic sim¬ 

ulations ( jManousakis & Walter|2015[ l. This provides some new 

insight on the inner stellar wind of HD 77581 and in particular 
on the f velocity law. The data and the simulations are described 
in Sect. [2]and[3j respectively. The comparison between data and 

simulation is presented in Sect. [3] and discussed in Sect. [4] Sec¬ 
tion 0 summarizes our results. 


2. Swift/BAT Data and Analysis 


The wide field of view of the Burst Alert Telescope (BAT, 
[Krimm et al.|2013| ) on board of the Swift satellite allows to mon¬ 
itor the complete sky at hard X-rays every few hours. For sources 
as bright as Vela X-l, lightcurves can be obtained in different en¬ 
ergy bands with a resolution of 1000 sec for almost a complete 
decade. 

The Swift/BAT reduction pipeline is described in |Tueller| 
et al. ( 2010| ) and Baumgartner et al. ( 2013| ). Our pipeline is based 
on the BAT analysis software HEASOFT v 6.13. A first analy¬ 
sis was performed with the task batsurvey to create sky im¬ 
ages in the 8 standard energy bands using an input catalogue of 
86 bright sources (that have the potential to be detected in sin¬ 
gle pointings) for image cleaning. Background images were then 
derived removing all detected excesses with the task batclean 
and weighted averaged on a daily basis. The variability of the 
background was then smoothed pixel-by-pixel using a polyno¬ 
mial model with an order equal to the number of months in the 
data set. The BAT image analysis was then run again using these 
averaged background maps. The image data were stored in a 
database organised by sky pixel (using a fixed sky pixel grid) 
by properly projecting the images on the sky pixels, preserving 
fluxes. This database can then be used to build images, spectra 
or lightcurves for any sky position. 

The result of our processing was compared to the standard 
results presented by the Swift team (lightcurves and spectra 
of bright sources from the Swift/BAT 70-months survey cata¬ 
logue 3 * * } and a very good agreement was found. 

The Swift/BAT lightcurves of Vela X-l were built in several 
energy bands. For each time bin and energy band a weighted mo¬ 
saic of the selected data is first produced and the source flux is 


1 http://swift.gsfc.nasa.gov/results/bs7Qmon/ 
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Fig. 1: The Lomb-Scargle periodogram based on the Swift/BAT 
14 - 100 keV energy band lightcurve. The red line shows the 
best derived orbital period, P = 8.964 ± 0.003 days. 

Table 1: Orbital periods and uncertainties. 


Time interval (MJD) 

Period (d) 

53370 - 56290 

8.964 ± 0.003 

53370 - 54830 

8.967 ± 0.007 

54830 - 56290 

8.965 ± 0.007 

53370 - 54344 

8.97 + 0.01 

54344-55315 

8.96 + 0.01 

55315 -56290 

8.97 + 0.01 


extracted assuming fixed source position and shape of the point 
spread function. The source signal to noise varies regularly be¬ 
cause of intrinsic variability, its position in the BAT field of view 
and distance to the Sun. The lightcurves span from ~ 53370 to ~ 
56290 MJD, i.e. over 300 orbital periods or almost 9 years. 

We have used the Lomb-Scargle (Press & Rybicki 1989) 
technique to determine the orbital period from the 14 - 100 keV 
Swift/BAT light-curve. The orbital period derived from the com¬ 
plete data-set is P 0 ,i=8.964 ± 0.003 days. Figure [I] shows the 
Lomb-Scargle power (black line), together with a Gaussian fit 
(in red). We searched for variations of the orbital period splitting 
the data in few equally spaced time intervals. The results are 
listed in Table [T] The orbital period remains constant within the 
uncertainties, and no secular evolution can be identified. The av¬ 
erage period together with the mid-eclipse reference time MJD 
= 53377.3964 have been used to obtain folded light-curves de¬ 
scribed below. 

The 14 - 100 keV lightcurve was also split in 8 successive 
periods (each spanning about a year) in order to study variabil¬ 
ity over these periods (Fig. |2j. Some major flaring activity can 
be identified in each of the periods. Vela X-l is known to be 
highly variable at hard X-rays and features major flares regu¬ 
larly ( |Kreykenbohm et ak||2008[ ). These flares (when averaged 
over a year) are not located at specific orbital phases, indicating 
that they are stochastic and much shorter than the orbital period. 

The averaged folded lightcurve shows a clear asymmetry be¬ 
fore and after the eclipse (Fig. [3]). This asymmetry is a likely 
signature of the accretion wake trailing the neutron star ( |Blondin| 
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Fig. 2: The folded lightcurve of Vela X-l in the 14 - 100 keV 
energy band, split in 8 consecutive datasets and shifted arbitrar¬ 
ily. The grey curves shown on the top and bottom include all 
available data for comparison. 


|et al.|!990, 1991 ). Several sgHMXBs ( |Manousakis et alj2012| 
and references therein) show a similar behaviour. The average 
absorbing column density of the accretion wake can be recon¬ 
structed from the ratio of the folded lightcurves observed on both 
sides of the eclipse. In Fig. [3] we compare the lightcurve preced¬ 
ing (</> = 0 - 0.5; black line) and following (</> = 0.5 - 1; red 
line) the eclipse. For the phases of <p ~ 0.5 - 0.8 the variations 
between the two curves are of the order of (1 - 2) x 10 23 crcr 2 
and related to source flares and to inhomogeneities in the stellar 
wind. A very strong flare, which occurred during the first year of 
observation, shows up close to phase 0.5. Shortly before ingress 
(i.e., <p ~ 0.8 - 0.9) the flux is significantly lower, corresponding 
to scattering in an additional absorbing column density of the 
order of 8 x 10 23 cm 2 . 

The BAT folded lightcurve during eclipse ingress and egress 
has an exceptional signal to noise. Besides constraining the ac¬ 
cretion wake, the lightcurve is particularly sensitive to the ve¬ 
locity of the unperturbed stellar wind close to the surface of the 
supergiant. An aspect that we can further probe comparing the 
observations to the results of simulations. 


3. Hydrodynamic Simulations 


The hydrodynamic simulations are largely based on our previous 
work (Manousakis & Walter|2015) but tuned to the new obser¬ 
vational results presented here. Only a brief description of the 
simulation code is provided here. The Euler equations dictate 
the motion of the wind, assuming mass, momentum, and energy 
conservation and account for the Roche potential and the line 
driven force ( |Blondin et al.|[1990l |1991| ).The mesh is fixed but 
non-uniform, allowing high resolution around the neutron star. 

The winds of hot massive stars are characterized observa- 
tionally by the wind terminal velocity and the mass-loss rate. 
Wind terminal velocities and mass-loss rates are typically in the 
range u m ~ 1000 - 3000 km s _1 and M w ~ Kr (6 ~ 7) M Q yr _I , 
respectively ( jKudritzki & Pulsj2000j i. The wind velocity profile 
is described by the ft- velocity law v = i/oo(l - R*/r)P, where /j is 
the velocity gradient (Castor et al .|1975) l. The values of Uco and 
[i are used to derive the CAR parameters characterizing the line 
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Fig. 3: Top panel. The observed folded light-curves preceding 
(black) and following (red) the eclipse. Bottom panel, inferred 
column density. The shaded vertical area indicates the eclipse. 
The light red horizontal area shows the artificial column density 
variability induced by the source short term intensity variability. 
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Fig. 4: The column density as a function of orbital phase derived 
from a smooth wind (dashed lines) and hydrodynamic simula¬ 
tions (solid lines) with a wind terminal velocity, ~ 1400 km 
s 1 (derived from 1-D simulation) and [3 — 0.5 (blue) and/3 = 0.8 
(green). 


driven force in our simulations using a 1-D simulation. These 
CAK parameters are then used in the 2-D code. 

The ionization of the wind by the neutron star is character¬ 
ized by the ionisation parameter l — Lx/nr? ls , where Lx is the 
average X-ray luminosity and n is the number density at the dis¬ 
tance r ns from the neutron star (Tarter et al. |1969| . Close to the 
neutron star, where > % crit ~ lCr 3 erg cm sec -1 , the ions re¬ 
sponsible for wind acceleration are highly ionized ( Kall man~&| 
McCray|1982|l and the radiative acceleration cuts-off. Although, 


the effects of the X-ray feedback on the wind acceleration force 
are complicated due to the large number of ions and line transi¬ 
tions contributing to the opacity (Abbott 1982; Stevens & Kail 


|man|1990] l a simplified approach allows the formation of a real¬ 
istic shock in front of the neutron star ( |Fransson & Fabian|1980| l. 

We have made use of VH- 1[^] code, described in depth by 
Blondin et al. (1990 |1991)>. The computational mesh employed 


in this study consists of 900 radial by 347 angular zones, ex 
tending from 1 to ~ 25 R* and in angle from -it to +n. The 
non-uniform grid is nested around the neutron star, reaching its 
highest resolution there, about ~ 10 9 cm . A co-rotating ref¬ 
erence system located around the center of mass is used. The 
outer boundary condition is characterized as an outflow while 
an absorbing boundary conditions (Hunt|197l Blondin & Pope 


|2009| is used at the stellar photosphere. The initial setup (i.e., 
wind density, velocity, pressure as well as the CAK parame¬ 
ters) is described in |Manousakis & Walter| ( |2015[ ) and results to 
/3 » 0.5 - 0.8 when the 2-D simulation stabilizes. 

We performed several simulation runs with a wind terminal 
velocity of » 1400, 1700 km s -1 ,/3 = 0.5, 0.8, and binary 
separations a = 1.76, 1.77, 1.78 R*. The mass loss rate of the 
donor star was always fixed to M « 4 x 10 -6 Mq. The time step 
of the simulations is ~ 1/10 sec. The simulations ran for roughly 
^ 3 orbital periods. The first ~ 3 days were excluded from the 
analysis, allowing the wind to reach a steady configuration. As¬ 
suming a mass to light conversion factor tj = 0.1, the resulting 
time-averaged mass accretion rate corresponds to an X-ray lumi¬ 
nosity in the range (2 - 8) x 10 36 erg/s. 


http://wonka.physics.ncsu.edu/pub/VH-l/ 


To compare the observations with the simulations we have 
built synthetic lightcurves derived for each simulation run. The 
neutron star is emitting as a point source and part of these X- 
ray photons are scattered in the stellar wind forming a diffuse 
source. The global intensity and geometry of the scattered emis¬ 
sion was calculated from the illumination and Thomson optical 
depth of each simulation cell. This is a reasonable assumption at 
hard X-rays. We assumed that the accretion wake extends verti¬ 
cally in the range |z|=0.5R* and that above these limits a smooth 
stellar wind applies. The very small flux (with a relative emis- 
sivity of 10 4 ) observed during eclipse originates from an ex¬ 
tended region of the wind and was accounted for. These synthetic 
lightcurves are shown as dashed lines in Fig. [5] 

We also calculated the column density Nh between the ob¬ 
server and the neutron star along the orbit cell by cell, excluding 
the highly ionised part of the wind, and corrected the synthetic 
lightcurve for the effect of scattering. These lightcurves are dis¬ 
played as continuous lines in Fig. [5] together with the Swift/BAT 
data (blue points). The correspondence between the simulations 
and the data, measured with a reduced^ 2 , are listed in Table [ 2 ] 
An excellent match is obtained for a steep (J3 = 0.5) velocity 
law and a terminal velocity of 1400 km/s. These values were the 
inputs to the 1-D simulation used to derive the CAK parameters 
and correspond to an unperturbed wind velocity of ~ 1400 km/s 
at the orbital radius and of 850 km/s at a distance of 0.2 R* of 
the stellar surface in the 2-D simulations. Binary separation sig¬ 
nificantly different from a — 1.77 R* induces large variations at 
eclipse ingress and egress that do not match the observed data 
Or >200 for 119 d.o.f.). 

The time-averaged absorbing column density derived from 
the simulations are shown in Fig. [4] for /3 = 0.8 (green) and 
/3 = 0.5 (blue) together with the predictions from an unperturbed 
smooth stellar wind (dashed lines). Significant Nh variability is 
expected from orbit to orbit. The average Swift/BAT spectrum 
for that phase range is not sensitive to low column densities so 
could not be used to further constrain the model parameters. 
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Fig. 5: The Swift/BAT 14 - 100 keV energy range folded light-curves of Vela X-l (blue points) together with the simulated eclipse 
profiles (dashed lines) using a binary separation a — 1.77 R*. Left: The light-curves corrected for scattering (solid line) were built 
assuming wind terminal velocity of Voo = 1400 and jj = 0.5. Right: The same as previous for the remaining models, vXXXXbYY, 
where XXXX is the wind terminal velocity in km s 1 and YY is the wind /? parameter. The residuals between the model and 
observations are shown in the bottom panel. They' 2 are listed in Table[2] 


Table 2: Simulation parameters (wind terminal velocity, beta pa¬ 
rameter and binary separation) and agreement (% 2 ) with the ob¬ 
served lightcurve. 


Model 

Voo (km s 2 ) 

P 

X 2 /d.o.f 

ar=1.76 R* 

vl708b®8 

1700 

0.8 

407/119 

vl488b®8 

1400 

0.8 

382/ 119 

vl780b®5 

1700 

0.5 

392/ 119 

vl480b®5 

1400 

0.5 

212/119 

a=1.77 R* 

vl780b®8 

1700 

0.8 

297/119 

vl48Qb®8 

1400 

0.8 

276/119 

vl780b®5 

1700 

0.5 

295/119 

vl480b®5 

1400 

0.5 

128/119 

£2=1.78 R* 

vl780b®8 

1700 

0.8 

452/119 

vl48Qb®8 

1400 

0.8 

405/119 

vl780b®5 

1700 

0.5 

399/ 119 

vl400b®5 

1400 

0.5 

253/ 119 


4. Discussion 

Spectroscopic observations of massive stars allow to constrain 
the mass loss rates and the terminal velocities of stellar winds 
( Puls et al.|2008[ and references therein). The P-Cygni line pro¬ 
files used for these estimates are formed far from the massive star 
at a few tens of stellar radii, where the /? parameter has a small 
effect. Orininally the models by Castor et al.|(|1975[> resulted in a 


steep ft = 0.5 law, while later on a shallower ft = 0.8 was intro¬ 
duced to match the observations ([Pauldrac h et al.|1986| |Friend| 
|& Abbott|1986[ |Puls et al.|2008[ ). lStevens ( |1993| > discussed the 
difficulties of low velocity gradients and the need for gradients 
varying with radius. In an attempt to study the stellar winds from 
Wolf-Rayet stars, Springmann (19941 employed a modified fi 
law mimicking a steep velocity field gradient (J3 = 0.5) close 
to the stellar surface and a smoother (J3 = 1.0) velocity gradient 
further away. Such a low /3 parameter close to the stellar surface 
was also found studying the soft X-ray absorbing column density 
along the orbit in 4U 1700-37 with EXOSAT (Haberl et al.ll989j > 
but these measurements were much more sensitive to ionisation 
than ours, based on hard X-ray data. 


The line driven instability predicts large density and velocity 
discontinuities which are not modelled in our simulations. Ac¬ 
cording to the latest models ( [Sundqvist & Qwocki|2013j >, this 
occurs at distances larger than one stellar radius from the stel¬ 
lar surface, i.e. well outside of the orbit in the case of Vela X-l, 
and should not play an important role in driving the X-ray vari¬ 
ability. Gravitation and the X-ray ionisation from the neutron 
star, that are taken into account in our simulations, remain the 
main drivers of the hydrodynamics. Photoionisation influences 
the wind within a few tens of the orbital radius (see also |Krtickaj 
|et al.|2012) but leaves the wind located at ±90° from the neutron 
star direction, responsible for the inress/egress scattering profile, 
largely unaffected. 


The hard X-ray scattering profile measured at eclipse ingress 
and egress is sensitive to the density of the stellar wind at a dis¬ 
tance of a few tenths of stellar radius from the stellar surface. At 
eclipse ingress it is in addition sensitive to the column density of 
the accretion wake formed close to the orbital radius. The den¬ 
sity of the wind depends directly on the /j parameter. A larger 
velocity gradient will generate a longer eclipse egress. Note that 
the wind velocity at 1.2 R , is about twice larger with fi = 0.5 
than with [i = 0.8 and has therefore a strong impact on the hard 
X-ray flux profile. The column density of the accretion wake also 
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Radius (RJ 


The simulation matches the observed data for a narrow set 
of parameters implying that the wind velocity close to the stellar 
surface is twice larger than usually assumed. The stellar wind in 
Vela X-l shows velocity field similar to these observed in some 
Wolf-Rayet stars (Springman n| 1994) . 

The hydrodynamic simulations of the wind of Vela X-l 
are still incomplete, however they provide already an excellent 
agreement with the observations for the smooth properties (this 
work) and the variability ( |Manousakis & Walter 20151. High- 
mass X-ray binaries provide informations on the initial acceler¬ 
ation of the stellar wind, which are unique and required to test 
models accurately. 
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5. Conclusion 

We have analysed almost a decade of S iv/'/f/BAT data of Vela X- 
1. Although variable, its hard X-ray emission, folded along the 
orbit, features a persistent asymmetric eclipse ingress and egress, 
measured with high signal-to-noise. This asymmetry originates 
from the scattering of hard X-rays around the neutron star at 
small and larger scales, including in the accretion wake. The X- 
ray variability is related to in-homogeneities in the stellar wind 
largely created by the gravity and ionisation of the neutron star 
itself ( [Manousakis & Walter|2015] >. 
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